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A new method of stochastically estimating the low-lying effects of quark propagation is proposed 
which allows accurate determinations of temporal correlations of single-hadron and multi-hadron 
operators in lattice QCD. The method is well suited for calculations in large volumes. Contributions 
involving quark propagation connecting hadron sink operators at the same final time can be handled 
in a straightforward manner, even for a large number of final time slices. The method exploits 
Laplacian Ifeaviside (LapH) smearing. Zn noise is introduced in a novel way, and variance reduction 
is achieved using judiciously-chosen noise dilution projectors. The method is tested using isoscalar 
mesons in the scalar, pseudoscalar, and vector channels, and using the two-pion system of total 
isospin 7 = 0, 1, 2 on large anisotropic 24'' x 128 lattices with spatial spacing as ~ 0.12 fm and 
temporal spacing at ~ 0.034 fm for pion masses « 390 and 240 MeV. 
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I. INTRODUCTION 

Recent discoveries of new hadronic resonances have 
generated much excitement in the field of hadron spec- 
troscopy. The current surge in experimental activity un- 
derlines the need for a better understanding of excited 
hadronic states from the theory of quantum chromody- 
namics (QCD). Presently, Markov-chain Monte Carlo es- 
timates of QCD path integrals defined on a space-time 
lattice offer the best way to make progress in this regard. 

Calculating the mass spectrum of excited-state hadron 
resonances is a key goal in lattice QCD. However, such 
calculations are very challenging. Computational limi- 
tations cause simulations to be done with quark masses 
that are unphysically large, leading to pion masses that 
are heavier than observed and introducing systematic er- 
rors in all other hadron energies. The use of carefully de- 
signed quantum field operators is crucial for accurate de- 
terminations of low-lying energies. To study a particular 
state of interest, the energies of all states lying below that 
state must first be extracted, and as the pion gets lighter 
in lattice QCD simulations, more and more multi-hadron 
states lie below the masses of the excited resonances. The 
evaluation of correlations involving multi-hadron opera- 
tors contains new challenges since not only must initial 
to final time quark propagation be included, but also 
final to final time quark propagation must be incorpo- 
rated. The masses and widths of resonances (unstable 
hadrons) cannot be calculated directly in finite- volume 
Monte Carlo computations, but must be deduced from 
the discrete spectrum of finite-volume stationary states 
for a range of box sizes [ll-Ill- 

Our approach to constructing hadron operators appro- 
priate for such calculations was outlined in Refs. [alll- 
Our first study of the nucleon and A excitations in the 



quenched approximation was presented in Ref. fj\, and 
nucleon results for two flavors of dynamical quarks ap- 
peared in Ref. [8]. A survey of excited-state energies 
in small volume for the isovector mesons and kaons us- 
ing Nf = 2 + 1 dynamical quarks was given in Ref. [l], 
along with results for the A, E, 5 baryons. Other recent 
progress in calculating excited-state energies in lattice 
QCD can be found in Refs. [iMl- AU of our results 
to date have been achieved in small volume with pi- 
ous having masses comparable to or heavier than about 
390 MeV. Our goal now is to extend our efforts into larger 
volumes and using lighter pions. To do this, the issue of 
multi-hadron states must be addressed. 

In this work, we focus on the problem of incorporat- 
ing multi-hadron operators into finite-volume excited- 
state spectrum calculations in lattice QCD. To compute 
the finite-volume stationary-state energies of QCD, one 
must first evaluate a matrix of temporal correlations 
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time-ordering, the source operators Oj(to) create the 
states of interest at an initial time to, and the sink opera- 
tors Oi{tp) annihilate the states of interest at a later time 
tp. The correlation functions Cij{t) can be expressed in 
terms of "path" integrals over quark ip, ip fields and gluon 
U fields involving the QCD action having the form 



(1) 



where K[U] is known as the Dirac matrix and S'g[C^] is 
the gauge-field action. Integration over the Grassmann- 
valued quark fields introduces a det K and factors of 
in the remaining integrals over the gluon U field, and 
when formulated on a Euclidean space-time lattice, such 
path integrals can be estimated using the Monte Carlo 
method with Markov-chain importance sampling. Incor- 
porating the det K in the Monte Carlo updating and eval- 
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uating the elements of (the quark propagators) are 
the most computationally demanding parts of the calcu- 
lations. 

Once estimates of a Hermitian matrix of tempo- 
ral correlation functions Cij{t) are obtained, sev- 
eral procedures for extracting the lowest stationary- 
state energies Eo, , E^,- ■ . in any given symme- 
try channel are available[16i [T7| . For example, let 
^n{t,to) denote the eigenvalues of the Hermitian matrix 
C{to)-^^^ C{t) C(^o)"^/^ where is some fixed reference 
time (typically small) and the eigenvalues, also known as 
the principal correlation functions, are ordered such that 
Ao > Ai > • • • as t becomes large. Then one can show 
that 

hm A„(t,io) - e-^"(*-*°). (2) 

Determinations of the principal correlators A„(i,to) for 
sufficiently large temporal separations t yield the desired 
energies 

The rows and columns of the gauge-covariant Dirac 
matrix K[U] can be viewed as compound indices which 
incorporate the lattice space-time site indices and the 
quark color, flavor, and spin indices. Hence, K is a very 
large matrix, and determining and storing all of the el- 
ements of is not possible. Symmetries are used to 
eliminate the need to compute all elements. Compu- 
tations are usually arranged such that the linear system 
of equations Kx = y needs to be solved for only a man- 
ageable number of source vectors y. For temporal cor- 
relations of single-hadron operators (excluding isoscalar 
mesons), invariance under all spatial and temporal trans- 
lations dramatically reduces the number of elements 
required. In such cases, the hadron creation operator 
needs to be considered only on one initial time slice and 
only at a single spatial site, yielding the so-called point- 
to-all method. A handful of points can be used to in- 
crease statistics. 

To study a particular eigenstate, the procedure by 
which energies are extracted from Monte Carlo estimates 
of temporal correlation functions using Eq. ^ requires 
that all eigenstates lying below the state of interest must 
first be extracted. As the pion gets lighter in lattice QCD 
simulations, more and more multi-hadron states will lie 
below the excited resonances, and multi-hadron opera- 
tors will be needed to accurately compute the energies 
of such states. For example, a good baryon-meson sink 
operator which annihilates a total zero momentum is typ- 
ically a superposition of terms having the form 

B{-p,t)M{p,t) = ^Y.^Bix,t)ipM{y,t)e^''-^''~y\ 

where V is the spatial volume of the lattice, 2p is the 
relative momentum, and tpBix^t) and (pMiy,t) are ap- 
propriate localized interpolating fields for a baryon and 
a meson, respectively. In the evaluation of the temporal 
correlations of such a multi-hadron operator, it is not pos- 
sible to completely remove all summations over the spa- 




FIG. 1: Examples of quark-line diagrams in multi-hadron 
correlators involving the i/) field on the later time t f connect- 
ing to a field also on the later time slice tp- The initial time 
is denoted by to. (Left) A two-meson correlator. (Right) The 
correlator of a baryon-meson system. 



tial sites on the source time slice using translation invari- 
ance. Hence, the need for estimates of the quark prop- 
agators from all spatial sites on a time slice to all 
spatial sites on another time slice cannot be sidestepped. 
Some correlators involve diagrams with quark lines orig- 
inating at the sink time tp and terminating at the same 
sink time tp (see Fig. [T|), so quark propagators involving 
a large number of quark-line starting times must also be 
handled. 

Finding better ways to incorporate the low-lying ef- 
fects of such slice-to-slice quark propagation for large 
numbers of quark source times is crucial to the success 
of our excited-state hadron spectrum project at lighter 
pion masses. A new method, known as distillation [l8|. 
uses a novel quark-field smearing procedure that facili- 
tates exact treatment of slice-to-slice quark propagation. 
Although distillation was found to work well, calculations 
with that method are costly, making it feasible only on 
small lattices. Here, we propose to combine the quark- 
field smearing used in Rcf. 18] with a new stochastic 
approach to estimating the quark propagators, resulting 
in a much more efficient treatment suitable for large vol- 
umes. Describing and testing this new method is the aim 
of this work. This method was briefly introduced with 
preliminary testing in Refs. [1, [l9l[20j. 

The remainder of this paper is organized as follows. 
The stochastic LapH method is described in Sec. [Hi 
Laplacian Heaviside quark-field smearing is reviewed, 
and our new stochastic treatment of quark propagation 
is detailed. The method involves Monte Carlo estima- 
tion of quark propagation using noise in the LapH 
subspace with variance reduction achieved through the 
introduction of suitable noise dilution projectors. The 
new method is compared to an earlier procedure which 
uses noise introduced on the space-time lattice itself. The 
number of inversions of the Dirac matrix needed in the 
new method is demonstrated to be insensitive to the vol- 
ume of the lattice. Details on how the temporal cor- 
relations of hadron operators are evaluated using the 
stochastic LapH method are then presented in Sec. IIIII 
Full source-sink factorization is seen to be another ad- 
vantageous feature of the method, especially for com- 
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putations of correlation matrices involving large sets of 
hadron operators. Various implementation details are 
given in Sec. IIVI Applications of the method to the 
isoscalar mesons in the scalar, pseudoscalar, and vec- 
tor channels and to the two-pion system of total isospin 
/ = 0, 1, 2 using anisotropic 24'^ x 128 lattices with pion 
masses m^r ~ 390 and 240 MeV are then presented in 
Sec. |Vl Conclusions are summarized in Sec. I VI I 



II. DESCRIPTION OF THE METHOD 

The use of smeared fields is crucial for successfully ex- 
tracting the spectrum of QCD in our Monte Carlo com- 
putations. Hadron operators constructed out of smeared 
fields have dramatically reduced mixings with the high- 
frequency modes of the theory that obscure the low-lying 
eigenstates of interest. Our operators are constructed us- 
ing spatially-smoothed link variables Uj (x) and spatially- 
smeared quark fields ipi^)- 

The spatial links are smeared using the stout-link pro- 
cedure described in Ref. [llj . Note that only spatial sta- 
ples are used in the link smoothening; no temporal sta- 
ples are used, and the temporal link variables are not 
smeared. 

The quark field for each quark flavor is smeared using 



(3) 



where x, y are lattice sites, a, b are color indices, a is a 
Dirac spin component, and the smearing kernel S is such 
that the smeared field behaves in exactly the same way 
as the original field under all time-independent symmetry 
transformations on a cubic lattice. For extracting ener- 
gies from temporal correlations, it is important that only 
spatial smearing is used. In other words, the smearing 
kernel is diagonal in time: Sab{x, y) oc Sx^y^- In addition, 
our smearing kernel is independent of spin. 

We use the new Laplacian Heaviside (LapH) quark- 
field smearing scheme which has been described in 
Ref. and is defined by 



S 



(4) 



where A is the three-dimensional gauge-covariant Lapla- 
cian defined in terms of the stout-smeared gauge field U, 
and CTs is the smearing cutoff parameter. The Laplacian 
matrix is given by 

3 

A'^^(x,2/;[/) = ^{c/f(a;)%,x + fc) 

k=l 

+UtiyrSiy,x- k) - 25{x,y)5'^''], (5) 

where x, y are lattice sites, and a, b are color indices. 
This is a Hermitian matrix which is block-diagonal in 
time. It is important to use the stout-smeared gauge 
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FIG. 2: The effective masses for temporal separation ts = 1 
for three representative nucleon operators against the LapH 
smearing cutoff crj. Results were obtained using Nf —2 + 1 
configurations on a 16^ x 128 anisotropic lattice with spacing 
as ~ 0.12 fm for stout-link smearing with = 10 iterations 
and staple weight ^ = 0.1. The circles show results (shifted 
downward by 0.04) for a single-site operator. The squares 
correspond to a singly-displaced nucleon operator, and the 
triangles are the results (shifted upward by 0.04) for a triply- 
displaced-T operator. The value al ~ 0.33 is observed to be 
a good choice. 



links when smearing the quark field since doing so dra- 
matically reduces the statistical errors in the correlations 
of the hadron operators we use which involve covariantly- 
displaced quark fields [8j. The gauge-covariant Laplacian 
operator is ideal for smearing the quark field since it is 
one of the simplest operators that locally averages the 
field in such a way that all relevant symmetry transfor- 
mation properties of the original field are preserved. 

Let Va denote the unitary matrix whose columns are 
the eigenvectors of A, and let Aa denote a diagonal ma- 
trix whose elements are the eigenvalues of A such that 



A = 1/a Aa Vl 



(6) 



The LapH smearing matrix is then given by 

S^VaQ{(tI+^a) Vl (7) 

All of the eigenvalues in Aa are real and less than zero. 
Hence, the matrix 0((t^ -I- Aa) has unit entries on those 
diagonal elements corresponding to eigenvalues whose 
magnitudes are less than a1 and zero entries for all other 
elements. Given that A is block-diagonal in time, each 
eigenvector has nonzero elements only on one time slice, 
so we can associate any given eigenpair with that particu- 
lar time. The eigenvalues of A always occur such that ap- 
proximately iVi, eigenvalues have magnitude smaller than 
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CT^ on each time slice. We have observed that the number 
of eigenvalues on each time slice that survive the Heav- 
iside function varies from time to time by only one or 
two in cases where Ny exceeds sixty or more. Hence, the 
Heaviside smearing matrix is well approximated by fix- 
ing Ny to the same value on each time slice and for each 
gauge configuration. 

Let Vs denote the matrix whose columns are in one-to- 
one correspondence with the eigenvectors associated with 
the Ny lowest-lying eigenvalues of —A on each time slice. 
Then our LapH smearing matrix is well approximated by 
the Hermitian matrix 



S = Vs vl 



(8) 



This is the actual smearing matrix used in our calcula- 
tions. Note that on a lattice having Nt time slices and 
Ns sites in each of the spatial directions, the matrix Vs 
has NyNt columns and NfN'^Nc rows, where Nc — 3 is 
the number of quark colors. The NyNt eigenvectors that 
form the smearing matrix span the so-called LapH sub- 
space. 

To set the parameter CTs, and hence Ny, several small 
simulations were done varying this parameter while com- 
puting the effective masses for a handful of simple meson 
and baryon operators. We chose the value of as that min- 
imized the effective masses at a chosen early time sepa- 
ration ts ■ The effective masses for ts — 1 for three repre- 
sentative nucleon operators are shown in Fig. [5] against 
values oi (7s- A single-site nucleon operator in which all 
three quark fields are taken at the same site is shown, 
as well as a singly-displaced nucleon operator in which 
one of the quarks is displaced away from the other, and a 
triply-displaced-T operator in which all three quarks are 
displaced from the others in a T configuration. The value 
as ~ 0.33 was chosen. This value is insensitive to which 
time interval is used as long as tg is small enough such 
that contributions from excited states have not decayed 
away. It is also insensitive to the choice of hadron oper- 
ator used and the quark mass. We expect as to change 
little with the lattice spacing. 

Evaluating the temporal correlations of our hadron op- 
erators requires combining matrix elements associated 
with various quark lines Q. Since we construct our 
hadron operators out of covariantly-displaced, smeared 
quark fields, each and every quark line in our computa- 
tion involves the following product of matrices: 

D'^^^Sn-^SD'-'''^^ (9) 

where fl — ^^K and D^'^ is a gauge-covariant displace- 
ment of type i. The displacement type can be trivial (no 
displacement), a displacement in a given single spatial di- 
rection on the lattice by some number of links (typically 
three) , or a combination of two or more spatial lattice di- 
rections. The use of $7 = liK is convenient for ensuring 
baryon correlation matrices that are Hermitian. 

An exact treatment of such a quark line is best accom- 
plished by writing 

Q = {V}Vl-~Ws) V}D^''^\ (10) 
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FIG. 3: The effect of the spatial lattice volume on the eigen- 
values of the gauge-covariant Laplacian operator. An is the 
n-th lowest eigenvalue of —A on a given time slice. The er- 
ror bars show the variation over different time slices and over 
a set of Nf =2 + 1 configurations. The lattice spacings Us 
are both near 0.12 fm, and the pion masses are both near 
0.70 GeV. Link smearing with = 10 iterations and staple 
weight ^ — 0.1 was used. For the 12'^ lattice, there are nine 
eigenvalues between 0.3 and 0.4, whereas for the 16^ lattice, 
there are 22 eigenvalues between 0.3 and 0.4, demonstrating 
that the density of eigenvalues is proportional to the spatial 
volume of the lattice (at sufficiently high values). The lowest- 
lying modes do not change very much with the lattice volume. 



then one needs to compute and store only the elements of 
the much smaller matrix V^Vl~^Vs instead of computing 
and storing a very large number of elements. Let 
N^ — 4 denote the number of Dirac spin components, 

and define yl!a°'\x) = T4(c, x; i) Sap, where a, /3 are spin 
indices, c indicates color, a; is a lattice site, and i refers 
to the column of Vs which is the i-th eigenvector of the 



Laplacian. Then solving the linear system Vlx = y 



(»,a) 

for X and all i,a by standard methods yields Vs ■ 
Hence, NyNfNd such inversions are required in order to 
obtain the full matrix V^n~'^Vs for each quark mass and 
each gauge configuration in the Monte Carlo ensemble. 
If only one source time slice is used in the hadron cor- 
relators, then NyNd inversions are required per quark 
mass per gauge configuration. Once multi-hadron oper- 
ators are included, however, sink-to-sink quark lines are 
needed, so NyNsnkNd inversions must be done, where 
Nsnk is the number of sink times. Generally, a hand- 
ful of hadron source times are used to improve statistics, 
so upon including multi-hadron operators, one finds that 
the number of inversions needed in practice ends up near 
NyNtNd. 

Solving the linear systems fix — y is a major compo- 
nent of the computational cost of evaluating the hadron 
correlators once a Monte Carlo ensemble is generated. 
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FIG. 4: The small effect of the light-quark mass on the eigen- 
values of the gauge- covariant Laplacian operator. A„ is the 
n-th lowest eigenvalue of —A on a given time slice. The error 
bars show the variation over different time slices and over a 
set of Nf — 2 + 1 configurations on a 16^ x 128 anisotropic lat- 
tice with Us ~ 0.12 fm for link smearing = 10 and ^ = 0.1. 



It turns out that the number of required eigenvec- 
tors on each time slice rises in direct proportion to the 
spatial volume of the lattice, as shown in Fig. [3] The 
number of eigenvectors is also fairly insensitive to the 
light quark mass, as shown in Fig. |4l Initial calculations 
on 16"^ lattices with spatial spacing as ~ 0.12 fm showed 
that Ny = 32 worked well. On 20^ lattices, Ny = 64 
was needed, and for the 24^ x 128 lattices used in this 
study, we found that Ny — 112 levels were below the 
CTg cutoff. We have generated gauge configurations on 
32'^ X 256 anisotropic lattices. On these lattices, we find 
that iV^, = 264, so the number of inversions needed be- 
comes NyNtNd > 270, 000 for each configuration and 
each quark mass, which is far too large to be feasible 
with current computing resources. 

Fortunately, an exact treatment of the quark lines is 
not needed. In fact, we have found that exact treatment 
of the quark lines is very wasteful. Given our use of the 
Monte Carlo method to evaluate the path integrals over 
the gauge link variables, the statistical errors in our es- 
timates of the hadron correlators are ultimately limited 
by the statistical fluctuations arising from the gauge-field 
sampling. Thus, we only need to estimate the quark lines 
to an accuracy comparable to the gauge noise from the 
Monte Carlo method. Such estimates can be obtained 
with far fewer inversions than required by an exact treat- 
ment of the quark lines. 

Random noise vectors 77 which satisfy E(r]i) — and 
^iViVj) = i^iii where E{) denotes an expected value as 
defined in probability theory, are useful for stochastically 
estimating the inverse of a large matrix as follows. As- 



sume that for each of Nr noise vectors, we can solve the 
following linear system of equations: ftX^' = rj^ for X^, 
where r labels the noise vectors r = 1, 2, • • • , Nr. Then 
X'' = n^^T]'', and E{Xi-q*) = so that a Monte Carlo 

estimate of O"-! is given by Vl'l « iV^' Ef=i ^z"??"- 
Unfortunately, this equation usually produces stochastic 
estimates with variances which are much too large to be 
useful. Variance reduction is done by diluting the noise 
vectors [23 - [2^ . A given dilution scheme can be viewed as 
the application of a complete set of projection operators 
P'^^y Define rj''V'\ = pC'^yy'', and define X' !''! as the solu- 
tion of = then a much better Monte Carlo 
estimate of il,^-^ is 

"^■^-]^EE^r''^?'*- (11) 

^ r=l b 

The dilution projections ensure exact zeros for many of 
the E{rjiri*) elements instead of estimates that are only 
statistically zero, resulting in a dramatic reduction in the 
variance of the 51^^ estimates. The use oi Zm noise en- 
sures zero variance in our estimates of the diagonal ele- 
ments E(r]irj*). The effectiveness of the variance reduc- 
tion depends on the projectors chosen. 

Earlier stochastic methods [2^, [2^ introduced noise in 
the full spin-color-space-time vector space, that is, on 
the entire lattice itself. However, since we intend to use 
Laplacian Heaviside quark-field smearing, an alternative 
is possible: noise vectors p can be introduced only in the 
LapH subspace. The noise vectors p now have spin, time, 
and Laplacian eigenmode number as their indices. Color 
and space indices get replaced by Laplacian eigenmode 
number. Again, each component of p is a random Zjv 
variable so that E{p) = and E{pp^) = Id, where Id 
is the identity matrix. Dilution projectors P^''^ are now 
matrices in the LapH subspace. In the stochastic LapH 
method, a quark line on a gauge configuration is evalu- 
ated as follows: 

= Eb-^(D^^'^5i7-iy,pWpp('=)v;pWp)t). (12) 

Displaced-smeared-diluted quark source and quark sink 
vectors can be defined by 

gV>\^p) = ^(^V.PC'V, (13) 

<^W(p) = D^^^Sn-^ VsP'^''^p, (14) 

and each quark line on a given gauge configuration can 
be estimated using 

, Nr 

Qir^ « if ^Ai.EE^l"(^^) ^^'^(pr, (15) 

^ r=l b 
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where the subscripts u, v are compound indices combin- 
ing space, time, color, spin, and quark displacement type, 
B is the flavor of the source field and A is the flavor of 
the sink field. The above quark line estimate has the 
form of an outer product expansion. Such estimates are 
frequently used in the compression of digital images, so 
the stochastic LapH estimate can be viewed as a lossy 
compression of the quark propagation. 

Occasionally, it is useful to estimate a quark line using 
75-Hermiticity to switch the source and sink. Using K'' = 
it is straightforward to see that another way to 
estimate a quark line is using 



Nr 



defining 



r=l b 



(16) 



(17) 



Eqs. (|T5|) and ([T6|) are meant to be used inside Monte 
Carlo estimates of path integrals over the gauge link vari- 
ables. To simplify matters, the Monte Carlo within a 
Monte Carlo computation can be combined into a single 
larger Monte Carlo calculation over both gauge link vari- 
ables and quark line noises, effectively setting Nji = 1 
for each gauge configuration. However, each quark line 
in a hadron correlator needs an independent noise to en- 
sure unbiased estimation. For example, a baryon correla- 
tor requires at least three noises per gauge configuration. 
Once inversions are done for a handful of such noise vec- 
tors for a given configuration, noise permutations can be 
used to increase statistics. 

The dilution projectors we use are products of time di- 
lution, spin dilution, and LapH eigenvector dilution pro- 
jectors. The full projector index b = {hT,bg,bi^) is a 
triplet of indices, where 6^ is the time projector index, 65 
is the spin projector index, and 6l is the LapH eigenvec- 
tor projector index. Our noise-dilution projectors have 
the form 



P. 



ib) 



tan: t'a'n' 



pibT) pibs) p{bL) 
t:t' a:a' n:n' 



(18) 



where t, t' refer to time slices, a, a' are Dirac spin in- 
dices, and n, n' are LapH eigenvector indices. For each 
type (time, spin, LapH eigenvector) of dilution, we stud- 
ied four different dilution schemes. Let N denote the 
dimension of the space of the dilution type of interest. 
For time dilution, N = Nt is the number of time slices 
on the lattice. For spin dilution, A'^ = 4 is the number 
of Dirac spin components. For LapH eigenvector dilu- 
tion, N — Ny is the number of eigenvectors retained on 
each time slice. The four schemes we studied are defined 
below: 



pib) _ 

p(b) „ X.. r . 
^ij — ^ij ^b, [Ji/Nj 



6 = 0, 
6 = 0, 
6 = 0, 



= ^ij h,i mod J 6 — 0, 



(no dilution) 
..,A^-1 (fuU dilution) 
, ., J - 1, (block- J) 
. . , J — 1, (interlace- J) 
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FIG. 5: Comparison of the new stochastic LapH method (tri- 
angles) with the earlier stochastic method using noise on the 
full lattice (squares) for the correlator C{t = 5) of a triply- 
displaced-T nucleon operator on a 16^ x 128 lattice. The 
vertical scale is the ratio of statistical error a (with no aver- 
aging over the six permutations of the three noises) over the 
error in the gauge- noise limit a^n, and in the horizontal scale, 
Nd is the number of Dirac-matrix inversions per source per 
quark line. Each point shows an error ratio using a particular 
dilution scheme. The LapH points lie significantly below the 
results from the lattice noise method, indicating a dramatic 
variance reduction. 



where i, j — 0, . . . , N — I, and we assume N/ J is an inte- 
ger. Note that each projector is a diagonal matrix with 
some or all of the diagonal elements set to unity and all 
other elements vanishing. We use a triplet (T, S, L) to 
specify a given dilution scheme, where "T" denote time, 
"S" denotes spin, and "L" denotes LapH eigenvector di- 
lution. The schemes are denoted by 1 for no dilution, F 
for full dilution, and B J and IJ for block- J and interlace- 
For example, full time and spin dilution 
LapH eigenvector dilution is denoted by 



J, respectively, 
with interlace-S 
(TF, SF, LIS). 

Introducing noise in this way produces correlation 
functions with significantly reduced variances, as shown 
in Fig. [5] Let C{t) denote the correlation function of 
a representative triply-displaced-T nucleon operator at 
temporal separation t. Let (Tgn represent the statisti- 
cal error in C(t = 5) using exactly-determined slice-to- 
slice quark propagators. In other words, (Tgn arises solely 
from the statistical fiuctuations in the gauge configura- 
tions themselves (the gauge noise limit). Let a denote the 
error in C{t — 5) using stochastic estimates of the quark 
propagators. The vertical axis in Fig. [5] is the ratio of 
the statistical error a in C{t = 5) over o-gn. Results are 
shown for a variety of different dilution schemes. In the 
lattice noise method, variance reduction is achieved with 
projectors which dilute in the time, spin, and color in- 
dices. Simple spatial dilutions are also used. The squares 
show results for dilution schemes with noise introduced 
in the larger spin-color-space-time vector space, and the 
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FIG. 6: Comparison of the new stochastic LapH method 
on 16'^ (triangles) and 20^ (squares) lattices for the corre- 
lator C{t — 5) of a triply-displaced-T nucleon operator on a 
16^ X 128 lattice. The vertical scale is the ratio of statistical 
error a (with averaging over the six permutations of the three 
noises) over the error in the gauge- noise limit o-gn, and in the 
horizontal scale, Nd is the number of Dirac-matrix inversions 
per source per quark line. Each point shows an error ratio 
using a particular dilution scheme. The number of Laplacian 
eigenvectors needed is 32 on the 16"^ lattice and 64 on the 
20^ lattice. The leftmost points correspond to the dilution 
scheme (TF, SF, LIS). For this scheme, a/a^n = 1.31 on the 
16^ lattice and a/agn ~ 1.32 on the 20'^ lattice. 



triangles show results for different dilution schemes using 
noise introduced only in the LapH subspace. The trian- 
gles show nearly an order of magnitude reduction in the 
statistical error, compared to the square symbols. The 
correlator for other time separations t and for a variety 
of other hadron operators were also examined. All of the 
observables we studied showed the same dramatic reduc- 
tion in the variance using the new LapH-noise method 
compared to the lattice-noise method. 

The number of Dirac matrix inversions needed in the 
stochastic LapH method to achieve a target statistical 
precision was found to be insensitive to the spatial vol- 
ume, despite the rapid increase in the number of LapH 
eigenvectors. Calculations on a 16^ and a 20^ lattice were 
carried out and the ratios cr/agn for various correlators at 
various time separations were compared. The error ratios 
for the representative triply-displaced-T nucleon correla- 
tor for time separation i = 5 on a 16^ lattice (triangles) 
are compared to those from a 20'^ lattice (squares) in 
Fig. E For the (TF, SF, LIS) dilution scheme, we found 
f/cgn — 1.31 for this quantity on the 16'^ lattice and 
f/cTgn = 1.32 on the 20^ lattice. Not only is the equality 
of these ratios on the two volumes remarkable, but also 
their closeness to unity is striking. Keep in mind that 
the number of Laplacian eigenvectors needed doubles in 
going from the smaller to the larger volume. These re- 
sults show that once a sufficient number of dilution pro- 
jectors are used, the number of inversions required by 



the stochastic LapH method does not increase with the 
lattice volume and are sufficient to essentially reach the 
gauge noise limit. Additional inversions of the Dirac ma- 
trix are totally unnecessary since they do not lower the 
error any further. Other time separations and a variety 
of other hadron operators were also studied and led to 
the same conclusions. 

Different dilution schemes were explored using 16"^, 
20'^, and 24"^ spatial lattices with spacing ~ 0.12 fm 
and light quark masses yielding pion masses ranging from 
240 MeV to 500 MeV, and we have found that the scheme 
(TF, SF, LIS) produces variances near that of the gauge 
noise limit for correlators which involve only quark lines 
that connect the source and sink time slices. Interlace- J 
and block- J were observed to work equally well for spin 
and LapH eigenvector dilutions. For correlators which in- 
volve quark lines that originate and terminate at the final 
sink time, the dilution scheme (TI16, SF, LIS) was found 
to work well. The interlacing in time enables us to evalu- 
ate quark lines that originate on any time slice. Results 
for several isoscalar correlators using (TI32, SF, LIS) on 
20 configurations were compared with (TI16, SF, LIS) 
and no differences in the variances were discernible, sug- 
gesting the gauge noise limit has essentially been reached. 

In the stochastic LapH method, the number of times 
that Qx = y must be solved is NpNp for each gauge- 
field configuration, where Np is the number of noises 
used and Np is the number of dilution projectors. Us- 
ing full time dilution (with four choices of source time 
to), full spin dilution and interlace-S LapH eigenvector 
dilution, then the to-to-tp (for sink time tp) quark lines 
require 128 inversions for each noise on each configura- 
tion. To accommodate a baryon-meson system, at least 
5 noises for these quark lines are needed. The tp-to- 
tp quark lines use interlace-16 time dilution, full spin 
dilution, and interlace-S LapH eigenvector dilution, re- 
quiring 512 inversions per noise per configuration. Two 
noises are required for these quark lines. Hence, a total 
of 5 X 12S -f 2 X 512 = 1664 inversions per configuration 
are needed to compute all baryons and mesons composed 
of u, d quarks. This number of inversions is the same 
for both the 24"^ and 32"^ lattices that we plan to use. 
An exact treatment of the quark propagation requires 
57,344 inversions per configuration on 24^^ x 12S lattices 
for = 112 and 270,336 inversions per configuration on 
32^ X 256 lattices for = 264. 



III. TEMPORAL CORRELATIONS OF 
HADRON OPERATORS 

Details on how the temporal correlations of hadron op- 
erators are evaluated using the stochastic LapH method 
are presented in this section. We limit our attention 
to four cases: baryon to baryon, meson to meson, two- 
meson to meson, and two-meson to two-meson systems. 
Other source to sink cases are straightforward general- 
izations of the four examples below. 
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A. Baryon to baryon correlations 

All of our hadrons are assemblages of basic building 
blocks which are covariantly-displaced, LapH-snieared 
quark fields: 



-A 

9aa, = V-aa 74 D 



(j)t 



(19) 



where a is a color index, a is a Dirac spin component, j is 
a displacement type, and ^ is a quark flavor. To simplify 
notation, the Dirac spin component and the displacement 
type are combined into a single index in what follows. 

Each baryon operator destroying a three-momentum 
p is a linear superposition of so-called elemental three- 
quark operators, which are gauge-invariant terms of the 
form 



abc q^^{x,t)q^p{x,t) q^^ {x,t). 



The "barred" three-quark elemental operators which cre- 
ate a momentum p have the form 



Sabc q%{x,t) q"p{x,t) qaa{x,t) 



We use hadron operators which transform irreducibly un- 
der all symmetries of the three-dimensional cubic lattice. 
Each baryon sink operator, being a linear superposition 
of the three-quark elemental operators, has the form 

Bi{t)^c%^'^tf^{t), (20) 

where I is a compound index comprised of a three- 
momentum p, an irreducible representation (irrep) A of 
the lattice symmetry group, the row A of the irrep, isospin 
and other flavor quantum numbers, and an identifier la- 
beling the different operators in each symmetry channel. 
The corresponding source operators are 

^^w = cS;$Xw- (21) 

The baryon correlation matrix elements are given by 

to 

where (. . .) denotes a vacuum expectation value, which is 
given by the usual ratio of path integrals over the fermion 
and gauge fields Wick-rotated into imaginary time. To 
simplify notation, we replace the average over all source 
times by a single fixed time tg, exploiting time-translation 
invariance, and obtain 



CaitF~to) = c%^^f_{ ^t^,^{t,) <p--^-ito) ). 

Expand the three-quark elemental operators in terms of 
the covariantly-displaced smeared quark fields. 



be 



'tp-[x — x) 



X ( qaaix, tp) qbp{x, tp) q^^{x, tp) 

X q^{x,ta) q§p{xM) q^{x,h 



where the three-momenta associated with I and I are as- 
sumed to be the same p, then evaluate the path integrals 
over the Grassmann fields to obtain a sum over products 
of quark lines, defining t = tp — to: 

XX 

QiAA)_Q(BB)_Q{CC)_ _ QiAA)_Q(BC)_Q{CB)_ 



QiAC)_g{BB)_g(CA)_ 

-^aa;c-i -^b0;bl3 '^cf,aa 



where time and spatial labels have been omitted, and 
(. . .)u is an expectation value defined as a ratio of path 
integrals over the gauge field U only, using the gauge field 
action and the fermion determinant as the path integral 
weight. Note that each quark propagator Q connects 
each source site x to each sink site x, as well as connect- 
ing color and spin components between the source and 
sink. Hence, the summations in the above equation are 
quite costly to carry out, and must be repeated over and 
over again for every pair of baryon operators. 

A dramatic simplification of the above equation can be 
achieved by using Eq. (1151) to estimate each quark line. 
The following quantity emerges as a key component of 
the resulting estimate: 



X<^LaL(Pl) ^titiP2) V^citiPzl (23) 

where 61,62,63 are noise dilution projector indices, and 
the short-hand notation ipk — ^{pk) has been used, where 
if is the quantity defined in Eq. (jl4[) . The baryon corre- 
lation matrix element is then given by 

Cfi{tp~t^) = (s^'^'^'l((pi,^2,^3;tF) 



ACE ^[616362] 



J ABC 



I 



{01,03, Q2;to) 



"ABC "7 (02, 01, 03, to) 



CBA ^[f>3f'2bi] 



{03, 02, 0i;to) 



"ABC 

"ABC \02, 03, 01, to) 



+ S^^^B\''^''^'^\g,,g„g,-,to)) ) (24) 

/ / U,p 

where S^§^ = Sad^be^cf and (. . .)u,p indicates an ex- 
pectation value over the gauge field U and any noises 
Pk ■ Again, the above equation uses the short-hand nota- 
tion ifk — ^{pk) and Qk — 0{pk), where the quark sinks 
Lp are defined in Eq. (1141) and the quark sources g are 
defined in Eq. (|T3l) . A,B,C are the quark flavors of the 
first, second, and third quarks as ordered in the B func- 
tions. To increase statistics, an average of the six permu- 
tations of the 1,2,3 superscripts labeling the quark lines 




FIG. 7: Graphical depiction of Eq. (|24|l for a baryon cor- 
relator C^jitp—ta) with source time to and later sink time 
tp- Each box represents a baryon function given by Eq. (|23p 
with the first quark located at the top of the box. Lines con- 
necting a Q with a </p indicate summation over their dilution 
projector identifiers. The same noise must be used at the two 
ends of any single line, and different noises should be used 
for different lines. Any line connecting quarks of different fla- 
vors represents a zero value. The asterisks indicate complex 
conjugation. 



can be used, and if the masses of all three quark lines are 
the same, this requires no further inversions of the Dirac 
matrix. 

A very useful feature of Eq. is the fact that the 
baryon correlator completely factorizes into a function 
associated vi^ith the sink time slice tp^ and another func- 
tion associated with the source time slice t^. Summations 
over color, spin, and spatial sites at the source have been 
completely separated from the color, spin, and spatial 
summations at the sink. The stochastic LapH method 
leads to complete factorization of hadron sources and 
sinks in temporal correlations, which greatly simplifies 
the logistics of evaluating correlation matrices involving 
large numbers of operators. Eq. ((24|) also shows that im- 
plementing the Wick contractions of the quark lines is 
also straightforward. Contributions from different Wick 
orderings within a class of quark-line diagrams differ only 
by permutations of the noises at either the source or the 
sink. In Eq. (1^^ , permutations of the noises at the source 
have been used since this is generally much less costly. 

Given the plethora of indices, a graphical representa- 
tion of the above equation is useful and is shown in Fig.[71 
The quark field tp is represented by a quark sink (/? or a 
p, and tp becomes & g ov & Tp. We represent a baryon 
given by Eq. p3|) by a box containing the quark sources 
or sinks vertically aligned with the first quark on the left 
in Eq. ([25]) located at the top of the box. We use lines 
connecting a g with a (or a 'g with a ^) to denote 
summation over the dilution indices associated with the 
connected g and Lp. The same noise must be used at the 
two ends of any single line, and different noises should be 
used for different lines. 



B. Meson to meson correlations 

Each meson operator destroying a three-momentum p 
is a linear superposition of quark-antiquark elemental op- 
erators which are linear superpositions of gauge-invariant 



(p- 
(p- 



"9 
-9 



9' 

1=0 



1=0 

:9_ 



"0 



"0 



FIG. 8: Graphical depiction of Eq. (I33p for a meson correlator 
C;j(tF— to) with source time to and later sink time tp- Each 
box represents a meson function given by Eq. (|32p with the 
flrst quark field located at the top of the box. Lines connecting 
a Q with a, ip 01 a 'g with a Tp indicate summation over their 
dilution projector identifiers. The same noise must be used 
at the two ends of any single line, and different noises should 
be used for different lines. Any line connecting quarks of 
different fiavors represents a zero value. The asterisks indicate 
complex conjugation. Contributions from the meson internal 
lines occur only for isoscalar mesons. 



terms of the form 

X 

(25) 

where q,q are defined in Eq. da,di3 are the spa- 

tial displacements of the q, q fields, respectively, from 
X, A,B indicate flavor, and a, /3 are compound indices 
incorporating both spin and quark-displacement types. 
The phase factor involving the quark-antiquark displace- 
ments is needed to ensure proper transformation prop- 
erties under G-parity for arbitrary displacement types. 
The "barred" operators which create a momentum p then 
take the form 

—AB 



W = ^e^P-(-+i('^=+''^))5,, q^^ix,t) qt{x,t) 



Each meson sink operator has the form 



(26) 



(27) 



(or is a flavor combination of the above form), where 
again, the / label includes the momentum p, the symme- 
try group irrep A, the row A of the irrep, and an identifier 
specifying the different operators in each symmetry chan- 
nel. The corresponding source operators are 



(28) 



The meson correlation matrix elements are given by 
CaitF-to) - ]^ E( ^^'(*^^) ^t(*o) )• (29) 

* to 

In terms of the elemental operators and using translation 
invariance, the above equation becomes, 

Caitp-to) = cSc|*( <i>^^{tp) $5(to) ), (30) 

using translation invariance to fix to a single to for no- 
tational convenience. Expand the elemental operators in 
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FIG. 9: Computation of the temporal correlation of a two- 
meson source at time to and a single-meson sink at time tp- 
The source mesons are assumed to be non-isoscalars. Each 
box represents a meson function given by Eq. (|32p with the 
first quark field located at the top of the box. Lines connecting 
a Q with a. (p OT a.'g with a Tp indicate summation over their 
dilution projector identifiers. The same noise must be used at 
the two ends of any single line, and different noises should be 
used for different lines. Any line connecting quarks of different 
fiavors represents a zero value. The asterisks indicate complex 
conjugation. The diagram with an internal line contributes 
only for isoscalar mesons. 

terms of the covariantly-displaced smeared quark fields: 

XX 

X q§^ix,to) q4six,to) ), 

where the three-momenta associated with / and I are as- 
sumed to be the same p. Next, the path integrals over 
the Grassmann fields are carried out, and one obtains, 
for t = tp — to, 

XX 

^gip-(^+i(d-+d^))/ p(£A) g{BB)_ 

omitting time and spatial labels. Eq. ^TE\\ or Eq. ^TE\\ 
can then be used to estimate the two quark propagators. 
In the first term, we find that it is advantageous to use 
Eq. ([T5| for the A quark line and Eq. ([T5|) for the B quark 
Unc. 

To proceed, define the following meson function: 

X 

^gtUpiT v^M), (32) 

where 61, 62 are noise dilution projector indices, and the 
short- hand notation (fk = '^{Pk) has been used again. 
The meson correlator is given by 

Cfi{tF~tQi) 




FIG. 10: Computation of the temporal correlation of a two- 
meson source at time to with a two- meson sink at time tp- 
All four mesons are assumed to be non-isoscalars. Each box 
represents a meson function given by Eq. p2|) with the first 
quark field located at the top of the box. Lines connecting 
a Q with a 95 or a ^ with a ^ indicate summation over their 
dilution projector identifiers. The same noise must be used at 
the two ends of any single line, and different noises should be 
used for different lines. Any line connecting quarks of different 
fiavors represents a zero value. The asterisks indicate complex 
conjugation. 

+ '5f|>^^'^'(^'l,v'l;^F)^<^''^k^2,&;^o)*) (33) 

where S'-^g = Sacked- The second term only contributes 
to isoscalar mesons. Again, color, spin, and spatial sum- 
mations at the sink have completely factorized from the 
summations at the source. This equation is graphically 
represented in Fig. HI 

C. More complicated correlations 

The graphical rules developed in the preceding sections 
can be applied to more complicated correlation matrix 
elements. The correlation of a two meson source with 
a single-meson sink is represented in Fig. [HI The source 
mesons are assumed to be non-isoscalars, otherwise there 
would be additional diagrams involving meson internal 
lines. The correlation of a two-meson source with a two- 
meson sink is illustrated in Fig. [TUl All four mesons are 
assumed to be non-isoscalars. We apply 75 Hermiticity 
only in cases where a V'(^o) at the source connects with 
a, tp{tp) at the sink. Full time dilution is the best choice 
for all quark lines connecting to and tp and to to to- For 
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all tp-to-tp quark lines, interlacing in source time must 
be used. 

To evaluate any correlation matrix element using the 
stochastic LapH method, one first must identify the vari- 
ous hadron functions that are needed and calculate them 
using Eqs. (1^ and ([5^ . These can be evaluated for a 
large set of hadron operators and stored on disk. The 
quark propagators are needed only at this stage of the 
computation. All color contractions and spatial sums 
are carried out in evaluating the hadron functions. Each 
hadron function for a given choice of noises takes up very 
little space on disk since each is an array over time and 
dilution indices only. The correlation matrix elements 
are then combinations of the different hadron functions 
for different noise selections. These final contractions in- 
volve only summations of dilutions indices. In this way, 
a large number of correlation matrices can be evaluated 
very efficiently. 



IV. IMPLEMENTATION DETAILS 

Our software is written in C-I--I- and links to the 
USQCD QDP-l-f /Chroma library [ill. Parts of our com- 
putations must be done using the full four-dimensional 
lattice, but other parts are best handled time slice by 
time slice in three dimensions. QDP-I--I- does not handle 
both three and four dimensional lattices simultaneously, 
so the different parts of the computations were done in 
separate runs using both 3d and 4d versions of our soft- 
ware. Special input/output routines were written to en- 
able 4d QDP++ to read and write 3d time slices of the 
lattice. 

Our computations are done as a sequence of tasks for 
each gauge configuration in the Monte Carlo ensemble. 
In the first task, the spatial links of the gauge configu- 
ration are smeared using the stout-link procedure. This 
task is done using a four-dimensional version of our soft- 
ware, but the smeared spatial links are written to disk 
as individual time slices suitable for input to the three- 
dimensional version of our software. In the second task, 
computation of the Laplacian eigenvectors is done time 
slice by time slice in three dimensions. In the third 
task, the eigenvectors for the different time slices are re- 
organized into four-dimensional eigenvectors correspond- 
ing to the different eigenvalues. The fourth task is the 
computation of the quark propagators. The inversions 
of the Dirac matrix must be done using the full four- 
dimensional lattice, but our results are written to disk 
once again as three-dimensional time slices. Formation 
of the hadron sources and sinks is accomplished in the 
fifth task using the three-dimensional version of our soft- 
ware. All of our hadron operators have definite three 
momenta which involve summations over all spatial sites 
of the lattice, so the resulting hadron sources and sinks 
are no longer lattice quantities. The final task is the 
assembly of the hadron sinks and sources to form the 
hadron correlation functions which can be accomplished 



using a serial version of our software. 

The eigenvectors of the gauge-covariant Laplacian are 
evaluated using a Krylov-Spectral Restarted Lanczos 
(KSRL) method which is a modification of the thick 
restarted Lanczos method described in Ref. [1^. Let A 
denote a Hermitian matrix whose lowest-lying or highest- 
lying eigenvectors are sought. Given a starting vector u, 
the KSRL method begins by constructing a Krylov space 
spanned by vectors u, Au, A^u, . . . , A^u. The submatrix 
of A defined in this basis is then diagonalized, and the 
eigenvalues and eigenvectors of this submatrix, known 
as the Ritz values and Ritz vectors, are approximations 
to those of the full matrix A. Convergence to the ex- 
act eigenpairs occurs as the Krylov space dimension in- 
creases, but a better procedure is to stop the growth of 
the Krylov space at some point, typically just above the 
number of desired eigenpairs, and restart the procedure 
using a different starting vector or vectors. The use of a 
certain number of Ritz vectors to restart the procedure is 
known as Krylov-Spectral restarting. Key issues in the 
method are determining how many Ritz vectors to use 
in restarting, determining the size of the Krylov space to 
use, and maintaining orthogonality of the Lanczos vec- 
tors in finite-precision mathematics. 

In our calculations, we use either a random vector or 
the vector whose components are all equal for the start- 
ing vector. Full global reorthogonalization is used at all 
steps. The decision to reorthogonalize multiple times is 
based on a simple criterion [29]: if the norm of the vector 
decreases by 1/k, where k = -s/S, then further reorthgo- 
nalization is done. A maximum of four reorthogonaliza- 
tions is enforced. Equation 5 in Ref. [30*1 is used to choose 
the number of Ritz vectors to keep, except that the num- 
ber must be at least as large as the number of converged 
vectors and cannot exceed the dimension of the Krylov 
space minus the number of converged and locked vectors 
minus twelve. For an approximate eigenvector x (with 
unit norm) and an estimate A of its corresponding eigen- 
value, the residual norm is defined by r = \\Ax — Xx\\. 
An eigenpair is considered converged when r < tol||^||, 
where tol is the desired tolerance and the matrix 2-norm 
is defined by ||A|| — max^j^o ||^2;||/||x||, and can be es- 
timated by the largest absolute value of any Ritz value 
encountered in the computation. 

In calculating the eigenvectors of A, Chebyshev accel- 
eration is used. The eigenvalues of —A are all real and 
lie between and some maximum value denoted by A^. 
We wish to determine the eigenvectors corresponding to 
the lowest-lying eigenvalues lying between and some 
cutoff Ac- The rate of convergence to solution increases 
with the spacing between the levels. Convergence is much 
faster for widely spaced levels. Hence, convergence can 
be accelerated by transforming the spectrum so that the 
desired part of the spectrum is more widely spaced. The 
following transformation is applied first: 

^-i+(Ar^(^+^-)- 
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FIG. 11: Masses mfit(t) obtained by fitting the correlators for single-site vr, Q., N operators to a cosh (exponential) form for the 
7r(r2, A'^) in the temporal range tmin to tmax. Results are shown for different tmin with tmax fixed to the value stated in the lower 
left corner of each plot. Open symbols indicate unacceptable fit qualities, and solid symbols show results with acceptable fit 
qualities Q. The top row shows results using 551 configurations of the 390 ensemble on a 24^ x 128 lattice, and the bottom row 
shows results with 584 configurations of the 240 ensemble on a 24^ X 128 lattice. The fit value given in each plot corresponds 
to the fit indicated by the red point. 



The above transformation maps the unwanted spectrum 
to the range — 1 • • • 1, and the desired part lies above 1. 
Chebyshev polynomials are then applied: 



A = T„(B). 



(35) 



Eigenvalues lying between -1 and 1 stay between -1..1, 
and the desired eigenvalues above 1 get spaced out to 
large and widely-separated values above 1. The lowest- 
lying eigenvalue of —A becomes the highest- lying eigen- 
value of A. Transforming the desired levels to the re- 
gion above 1 is most convenient since it allows the use 
of Chebyshev polynomials of any order, both even and 
odd. The Chebyshev polynomials are applied using the 
following recurrence relation: 



Tn{x) 



1, 



Ti{x) 



(36) 
(37) 



For calculations done on our anisotropic 24'^ x 128 lat- 
tices, we need to compute the lowest-lying Ny = 112 
eigenvectors on each time slice. A Krylov space dimen- 
sion of 160 was found to work well, and Xl — 15 and 



Ac = 0.5 were appropriate. Chebyshev polynomials of 
order 8 were used, and the residual tolerance was set to 
10~^. Convergence of all iV„ levels occurred within a 
dozen or less restarts. 

The LapH eigenvectors are uniquely determined only 
to within an overall phase. Given the way in which Zn 
noise is injected in the LapH subspace, one sees that 
a given quark line is not invariant under a change of 
the phase multiplying each eigenvector (due to the off- 
diagonal pieces not being exactly zero). It turns out that 
changing the phase is equivalent to changing the noise by 
a U(I) phase. This is not a problem, but erroneous re- 
sults can occur if the original eigenvector files used to de- 
termine the quark sinks get deleted and the eigenvectors 
have to be reconstructed for making the hadrons. With 
different run parameters, the eigensolver could produce a 
different phase. The introduction of a phase convention 
eliminates this potential problem. 

Once the needed eigenvectors of the Laplacian are com- 
puted and stored, the next step is to compute all ele- 
ments of 17" ly^pWp'' ■ There are only NtNy elements 
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FIG. 12; Products of at and the nucleon and fi-baryon masses 
against {171,, /mnf' for fixed fi = 1.5, rUs = —0.0743 and vary- 
ing rrin = rud- The two leftmost points for each baryon 
are from this work, and the three rightmost points are from 
Ref . [31I ] . The vertical dashed line indicates the physical value 
of [m-n/maY- 



to store for each noise r and each dilution projector &, so 
storage of these quark propagation coefficients is modest. 
Disk storage is actually dominated by the LapH eigen- 
vectors. Another nice feature is the fact that the quark 
propagation coefficients are gauge invariant, as long as 
the eigenvector phases are handled appropriately. Solv- 
ing fix = y for X with y = VsP^^'^p^ is accomplished 
using a mixed-precision improved version of the bicon- 
jugate gradient method with even-odd preconditioning. 
This was found to be the fastest inverter available in 
Chroma. Occasionally convergence is not achieved, and 
a slower conjugate gradient solver is applied to the sys- 
tem ^l^^x = ri'l^y. 

Our correlator estimates and their variances are insen- 
sitive to the value of N used for the noise, as long 
as N is not too small. We found that = 4 produced 
results indistinguishable in quality from those of larger 
N . Hence, we use Z4 noise in all of our computations. 
We identify a Z4 noise vector for an ensemble of gauge 
configurations by a 16-bit unsigned integer s. To create 
a noise vector p'*' for a gauge configuration labeled by 
an RHMC trajectory number k (assumed to have a value 
ranging from to 2^^ — 1), a 32-bit unsigned integer ra is 
first formed in a particular manner using the 16 binary 
digits of s and the 16 bits of k. Although the procedure of 
forming m is arbitrary, the same procedure must be used 
in every instance. The 32-bit unsigned integer m is then 
taken as a seed to the 32-bit Mersenne twister random 
number generator which is used to create the Z4 noise 



i,a) for each LapH eigenvector, labeled by time 
t and level i, and for each spin index a. The elements 
of p^^^ are generated in a particular order that is always 
the same. Each Z4 element is chosen using the sequence 
of bits obtained from the current state of the Mersenne 
twister, taking two bits at a time. It was found that the 
linear congruential generator in QDP-|--f /Chroma is not 
adequate for generating the Z4 noise and leads to serious 
errors in some instances. 



V. INITIAL APPLICATIONS 

Our initial development of the stochastic LapH method 
was done using a small 16'^ spatial lattice which is not 
very interesting for hadron physics. Since the main rea- 
son for pursuing the method is to apply it on large lattices 
for both single-hadron and multi-hadron correlators, we 
proceeded to test the method by studying several simple 
hadronic systems requiring sink-to-sink quark lines on a 
reasonably large 24'^ x 128 anisotropic lattice having spa- 
tial volume (3 fm)^. 

Two ensembles of gauge configurations were used. 
These ensembles were generated using the Rational 
Hybrid Monte Carlo (RHMC) algorithm (31), which is 
a variant of the hybrid molecular-dynamics (HMC) 
algorithm [ss'l suitable for Nf = 2 -|- 1 quark fiavors. The 
updating algorithm is a Metropolis method with a so- 
phisticated means of proposing a global change to the 
gauge and pseudofermion fields. A fictitious momentum 
is introduced for each link variable with a Gaussian dis- 
tribution and Hamilton's equations involving these mo- 
menta and the original action as a potential energy are 
approximately solved for some length of fictitious time, 
known as an RHMC trajectory. An improved anisotropic 
clover fermion action and an improved gauge field ac- 
tion were used[3l[. In both ensembles, /? — 1.5 is used 
and the s quark mass parameter is set to nis — —0.0743 
in order to reproduce a specific combination of hadron 
masses [Slj. In one ensemble, the light quark mass pa- 
rameters are set to niu = rrid = —0.0840 so that the pion 
mass is around 390 MeV using one particular way of set- 
ting the scale, discussed below. In the other ensemble, 
niu = trid — —0.0860 are used, resulting in a pion mass 
around 240 MeV. We refer to these ensembles as the 390 
and 240 ensembles, respectively. 

We calculated the masses of the pion, the nucleon, and 
the ri-baryon. Our results are shown in Fig. [TT] for the 
two ensembles on a 24^^ x 128 lattice. This figure demon- 
strates that the use of our stochastic estimates of the 
smeared quark propagators still leads to high accuracy 
results of standard quantities. The nucleon and fi-baryon 
masses times at are shown in Fig. 1121 against (m-n/ran)^- 
Results from Ref. [3l| are also included in this figure. 
Our goal in this work is simply to present and test the 
stochastic LapH method, so we defer a detailed analy- 
sis of these results until a later publication. However, 
it is encouraging that fitting the three leftmost fi-baryon 
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FIG. 13: Correlators C{t) against temporal separation t for single-site operators which produce the isoscalar pseudoscalar 
(PS), vector (V), and scalar (S) mesons. Results in the top row were obtained using 210 configurations (135 for the scalar 
channel) of the 390 ensemble. Results in the bottom row were obtained using 198 configurations of the 240 ensemble. In the 
legends, "fwd" refers to contributions from the diagram containing only forward-time source-to-sink quark lines, "smt" refers 
to contributions from the diagram containing only quark lines that originate and terminate at the same time. For the scalar 
channel, the "smt" contribution has a vacuum expectation value subtraction. Forward-time quark lines use dilution scheme 
(TF, SF, LIS) and same-time quark lines use (TI16, SF, LIS). The lattice size is 24^ x 12S for all the results shown here. 



points to a form linear in (tott/^o)^ and fitting the three 
leftmost nucleon points to an empirical form linear in 
'niTr/mn yields mN/mo, « 0.556 at the physical value of 
rriTj/mn, which compares well with the observed 0.561 
value. 

Our calculations determine all hadron masses in terms 
of the temporal lattice spacing a* . In order to express the 
hadron masses in terms of MeV, a value of a^^ must be 
specified using an appropriate renormalization scheme. 
Away from the physical point, different renormalization 
schemes will lead to different choices of a^'^. One partic- 
ular scheme that has been used in the past uses the mass 
of the rj baryon to set the scale when the strange quark 
mass is close to the value that reproduces the physical 
value of (2m|- — m^)/mQ. Using this scheme, we find 
an inverse temporal spacing a^^ = 5.661(17) GeV for 
the 390 ensemble and a^^ = 6.015(17) GeV for the 240 
ensemble. Since the ratio of spatial spacing over tempo- 
ral spacing has been tuned to a value near 3.5, we have 
Us ~ 0.12 fm for both of these ensembles. Our values for 



the pion and nucleon masses are to^ = 0.3911(14) GeV 
and TOjv = 1.1781(58) GeV on the 390 ensemble, and 
= 0.2439(20) GeV and m^v = 1.048(14) GeV on the 
240 ensemble. An alternative scale-setting scheme would 
be to extrapolate the r2-baryon mass results for different 
TO„,md but fixed /3,ms to the physical value of m-^jmi-i 
using a form linear in (jnT^/m^)'^ motivated by heavy 
baryon chiral perturbation theory, then use the VL mass 



to determine 



Doing this yields ~ 6.3 GeV and 



pion masses 250 and 430 MeV for our two ensembles. 

Results for the isoscalar mesons in the pseudoscalar, 
vector, and scalar channels and the two-pion system of 
total isospin 7 = 0,1,2 are presented in Figs. [131 HH HH 
and [ini In these results, the dilution scheme (TF, SF, 
LI8) is used for all quark lines connecting source time 
to the sink time tp and to to t^. Four widely-separated 
source times to were used on each gauge configuration. 
For all tp-to-tF quark lines, the dilution scheme (TI16, 
SF, LI8) is used. Observables are evaluated using con- 
figurations separated by ngcp RHMC trajectories, where 
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FIG. 14: Masses mflt(t) obtained by fitting the correlators C{t) shown in Fig. [T3]to a cosh form in the temporal range t-mi-n to 
tmax. Results are shown for different tm\-n with tmax fixed to the value stated in the lower left corner of each plot. Open symbols 
indicate unacceptable fit qualities, and solid symbols show results with acceptable fit qualities Q. The top row corresponds to 
the 390 ensemble, and the bottom row corresponds to the 240 ensemble. The left-hand-side plots show results for the 77 and 
TT pseudoscalar mesons, and the right-hand-side plots show results for the uj and p vector mesons. The scalar channel is not 
shown here since a reliable extraction of the lowest-lying energy in this channel needs a two-pion operator. The lattice size is 
24^ X 128 for all the results shown here. 



TT-scp = 20 for the two-pion correlators and rigcp = 40 for 
the isoscalar meson correlators. Jackknife binning shows 
autocorrelations to be suitably small. 

Our goal here is simply to test the stochastic LapH 
method, so simple single-site operators involving only the 
light u, d quarks are used for the isoscalar mesons, and 
single-site pion operators are used to make the two-pion 
states with zero and non-zero relative momenta. The 
temporal correlations of such simple operators have sig- 
nificant contaminations from higher-lying states, so that 
the effective masses associated with these correlations 
tend to a plateau rather slowly. Future work will make 
use of more sophisticated spatially-extended operators. 
The issue of mixing with Is operators is not addressed 
in these tests, and no vacuum-expectation- value subtrac- 
tion is used for the 77 correlator. In chirally-symmetric 
fermion formulations, the expectation value of the un- 



smeared, isosinglet pseudoscalar operator would be pro- 
portional to the topological charge, which has notoriously 
long autocorrelation times and may not be sampled prop- 
erly in a Monte Carlo simulation [s^. fssj . This can show 
up as a non-zero vacuum expectation value for the 77, 
which disappears as the volume increases. Our test re- 
sults do not take such effects into account, but future 
work will investigate this. 

In Fig. [131 the contributions to the isoscalar tempo- 
ral correlations C{t) from the diagram containing only 
forward-time source-to-sink quark lines are shown with 
label "fwd" , and the contributions from the diagram con- 
taining only quark lines that originate and terminate at 
the same time are shown with label "smt". The total 
correlator is also shown in each case. In the vector chan- 
nel, the contribution from the same-time diagram is very 
small and the total correlator can barely be distinguished 
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FIG. 15: Correlators C{t) against temporal separation t for two-pion operators with total isospin 7 = 0, 1, 2 and zero total 
momentum. S'-wave results have zero relative momentum, P-wave has minimal non-zero on-axis relative momenta. Results were 
obtained using 584 configurations of the 240 ensemble. In the legends, "fwd" refers to contributions from diagrams containing 
only forward-time source-to-sink quark lines, "smt" refers to contributions from diagrams containing only quark lines that 
originate and terminate at the same time, and "box" refers to diagrams containing both kinds of quark lines. Forward-time 
quark lines use dilution scheme (TF, SF, LIS) and same-time quark lines use (TI16, SF, LIS). The lattice size is 24'^ x 128 for 
all the results shown here. 



from the forward-line diagram contribution, so the "fwd" 
contribution is not shown. In the scalar channel, the ac- 
curacy of the "smt" contribution is particularly remark- 
able since a large vacuum expectation value has been 
subtracted. 
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FIG. 16: Energies Efit{t) obtained by fitting the correlators 
C(t) shown in Fig. [15] to a cosh 4- constant form in the tem- 
poral range tmin to fmax. Results are shown for different t-mi-n 
with fmax fixed to the value stated in the lower left corner of 
the figure. Open symbols indicate unacceptable fit qualities, 
and solid symbols show results with acceptable fit qualities 
Q. These results were obtained using 5S4 configurations of 
the 240 ensemble. The horizontal dashed lines indicate the 
energy of two free pious at rest. The lattice size is 24'^ x 128. 



The correlators in Fig. [T3l were used to extract various 
isoscalar meson masses. The pion and p masses can be 
obtained from the forward-line contributions to the pseu- 
doscalar and vector correlators, respectively. Correlated- 
fits to A(e~™* -|- e~'"(''^*~*)) for temporal separations 
between imin and tmax were done to extract the masses 
of the particles. Jackknife sampling was used to estimate 
the data covariance matrix, and bootstrap sampling was 
used to compute the uncertainties in the fit parameters. 
Results are shown in Fig. [T3] for various imin values, with 
imax fixed to the value stated in each plot. Even us- 
ing such simple hadron operators, fairly accurate mass 
extractions are obtained. Future use of better opera- 
tors will certainly improve these results. Results are not 
shown for the scalar channel since the lowest-lying energy 
in this channel is a two-pion state. Extractions of the 
energies in the scalar channel are best done with a corre- 
lator matrix using both single-hadron and two-pion oper- 
ators. The excellent statistical precision obtained for the 
correlators at small temporal separations suggests that 
diagonalizations of future correlation matrices estimated 
with stochastic LapH will be stable and accurate. 

With current Monte Carlo algorithms on presently 
available computing resources, it remains impractical to 
use light u, d quark masses tuned to properly reproduce 
the pion mass. Hence, the u, d quark masses used here 
yield a pion mass which is too heavy, making compar- 
ison to experiment somewhat problematical. Using the 
ri-baryon mass to set the inverse temporal spacing, we 
find masses = 576(59) MeV, rup = 820(13) MeV, 
and m„ = 863(21) using 198 configurations of the 240 
ensemble. The experimental values are = 548 MeV, 
TOp = 776 MeV, and to„ = 783 MeV. Future work will 
use better operators and all 584 configurations to achieve 
improved results. 
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Our results for the energies of two light pions are shown 
in Figs. [T5]and[T51 We studied S'-wave states of zero rel- 
ative momentum and total isospin 7 = and 7 = 2, 
as well as a P-wave with minimal non-zero on-axis rela- 
tive momenta in the I — 1 channel. In Fig. I15[ contri- 
butions to the correlators from the diagrams containing 
only forward-time source-to-sink quark lines are labeled 
by "fwd" , contributions from diagrams containing only 
quark lines that originate and terminate at the same time 
are shown as "smt" , and contributions labeled by "box" 
are those from the diagrams containing both kinds of 
quark lines (see Fig. [TU]). Energies were extracted using 
correlated-x^ fits to the form A + B(e~^* + e-^^^'"*)) 
in the range tmin to tmax- The results for different tmin 
are shown in Fig. [161 for ^max fixed to the value stated 
in the figure. Open symbols indicate unacceptable fit 
qualities, whereas solid symbols indicate results from fits 
of acceptable quality Q. The constant term in the fit 
form arises from a source pion propagating forward in 
time interacting with a sink pion propagating backwards 
in time and other similar contributions. The constant 
term was clearly evident in the 1 = 2 channel, but was 
consistent with zero in the / = channel. Hence, the 
1 = results shown in Fig. [15] were done setting the con- 
stant term to zero. This figure demonstrates that the 
stochastic LapH method can provide sufficient accuracy 
to see the difference of these two-pion energies from the 
energy of two free pions at rest, indicated by the horizon- 
tal dashed lines. In the 7 = 1 channel, the p-meson is ex- 
pected to be the lowest-lying energy level, so a correlator 
matrix including single-hadron and two-pion operators is 
necessary to reliably extract the low-lying spectrum in 
this channel. This will be done in future work. Again, 
we emphasize that only very simple operators were used 
here, and future use of better operators will improve the 
accuracy of these energies. 

VI. CONCLUSION 

A new method of stochastically estimating the low- 
lying effects of quark propagation was proposed which 
allows accurate determinations of temporal correlations 
of single-hadron and multi-hadron operators in lattice 
QCD. The method enables accurate treatment of hadron 
correlators involving quark propagation from all spatial 
sites on one time slice to all spatial sites on another time 
slice. Contributions involving quark lines originating at 
the sink time tp and terminating at the same sink time tp 
are easily handled, even for a large number oi tp times. 



The effectiveness of the method can be traced to two of 
its key features: the use of noise dilution projectors that 
interlace in time and the use of Zm noise in the subspace 
defined by the Laplacian Heaviside quark-field smearing. 
Introducing noise in the LapH subspace results in greatly 
reduced variances in temporal correlations compared to 
methods that introduce noise on the entire lattice. Al- 
though the number of Laplacian eigenvectors needed to 
span the LapH subspace rises dramatically with the spa- 
tial volume, we found that the number of inversions of the 
Dirac matrix needed for a target accuracy was remark- 
ably insensitive to the lattice volume, once a sufficient 
number of dilution projectors were introduced. 

In addition to increased efficiency, the stochastic LapH 
method has other advantages. The method leads to com- 
plete factorization of hadron sources and sinks in tempo- 
ral correlations, which greatly simplifies the logistics of 
evaluating correlation matrices involving large numbers 
of operators. Implementing the Wick contractions of the 
quark lines is also straightforward. Contributions from 
different Wick orderings within a class of quark-line di- 
agrams differ only by permutations of the noises at the 
source. 

The method was tested using the isoscalar mesons in 
the scalar, pseudoscalar, and vector channels, and using 
the two-pion system of total isospin / = 0,1,2 on large 
anisotropic 24'^ x 128 lattices with pion masses m^r ~ 390 
and 240 MeV. Given the success of these tests, we are 
now applying the stochastic LapH method to compute 
the excitation spectrum of both mesonic and baryonic 
stationary-states of QCD in large finite volume. 
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